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Abstract. We present a simple and easy to apply methodology for using high-level self-submitting 
parallel job queues in an MPI environment. Using CH — h, we implemented a library of functions, 
MPQueue, both for testing our concepts and for use in real applications. In particular, we have 
applied our ideas toward solving computational combinatorics problems and for finding bifurcation 
diagrams of solutions of partial differential equations (PDE). Our method is general and can be 
applied in many situations without a lot of programming effort. The key idea is that workers 
themselves can easily submit new jobs to the currently running job queue. Our applications involve 
complicated data structures, so we employ serialization to allow data to be effortlessly passed 
between nodes. Using our library, one can solve large problems in parallel without being an expert 
in MPI. We demonstrate our methodology and the features of the library with several example 
programs, and give some results from our current PDE research. We show that our techniques are 
efficient and effective via overhead and scaling experiments. 
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Our methodology is mainly motivated by our work in computational combinatorics |27| and 
partial differential equations (PDE) \22 \ [23 l [2~i]. In combinatorial applications, we have the needs 
and methods to follow trees. The bifurcation diagrams we study in PDE are similar in structure 
to these trees, and so we can supervise their generation using similar ideas. Traditional approaches 
to parallel PDE solvers obtain speed up by distributing operations on large matrices across nodes, 
"^j- , Our methodology instead uses job queues to solve such problems at a higher level. Our parallel 

^1 ' implementation for creating bifurcation diagrams for PDE is suggested by the nature of a bifurcation 

diagram itself. Such a diagram contains branches and bifurcation points. These objects spawn each 
other in a tree hierarchy, so we treat the computation of these objects as stand-alone jobs. For 
PDE, this is a new and effective approach. Processing a single job for one of these applications 
may create several new jobs, thus a worker node needs to be able to submit new jobs to the boss 
^ | node. These jobs require and produce a lot of complicated data. For our proof of concept tests 

and real applications, we have implemented our ideas in a library of C++ functions, in part to take 
advantage of that language's facility with complicated data structures. In this article, we refer to 
our library as MPQueue. 

The Message Passing Interface (MPI) [12] is a well-known portable and efficient communication 
protocol which has seen extensive use. In spite of its advantages, MPI is difficult to use because 
it does not have a high-level functionality. MPQueue is small, lightweight, and remedies some 
of these shortcomings. It uses MPI behind the scenes, is easy to use, and is Standard Template 
Library (STL) friendly. MPQueue has a high-level functionality allowing for the rapid development 
of parallel code using our self-submitting job queue technique. In order to easily pass our data 
structures between nodes, our library uses the Boost Serialization Library (BSL) |25| . It was a 
design goal of ours to make everything as simple as possible, so that one does not have to worry 
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about the details of MPI or serialization when applying our methods. MPQueue is freely available 
from this paper's companion website |21| , 

Section [2] contains a brief overview of some popular software packages that facilitate parallel 
programming. These fall into roughly three categories, threaded applications for shared memory 
multi-core processing, message passing for multi-node, distributed memory clusters, and single- 
system image clusters. Our library uses message passing and is intended for use on distributed 
memory clusters. Sections [3] and 0] explain line- by- line two simple example programs, one for 
factoring an integer, and the other for squaring a matrix. The purpose of these examples is not to 
provide efficient solutions for these demonstration applications, but to present our methodology and 
describe the key features of MPQueue in full detail. In Section [5] we give a more serious example, 
namely that of finding placements of non-attacking queens on a chessboard. This code is more 
complicated, and while not optimal, it does solve the problem on a 20 x 20 board, which would not 
be possible in a timely manner using serial code. We also investigate efficiency, scaling, and speedup 
for this example. A computationally intensive and mathematically complex example result for PDE 
can be found in Section El The full mathematical details of the particular PDE we are interested 
in can be found in [24] . along with our state of the art numerical results which use MPQueue. 
In Section we give an overview of the implementation of MPQueue. This section also includes 
an overhead experiment and summarizes the evidence of our implementation's solid performance. 
Section [8] contains some concluding remarks, including possible future refinements to MPQueue. 

The companion website |21| for this paper contains: the MPQueue source code, the library 
reference guide (see also Table IT. 1 P , all the files needed to compile and execute the three example 
programs, and for the less initiated, line-by-line descriptions of the three example programs. 

2. Related work 

In this section we give a short description of some existing systems offering parallel computing 
solutions. 

2.1. Shared memory multithreading systems on multicore processors. Many existing li- 
braries share some features with MPQueue, but unlike MPQueue, use a shared memory model 
with multithreading. These systems are not designed for use with a distributed memory cluster. 
They could conceivably be effectively used with some single-system image cluster software (see 
Section |2~3|) . 

• Intel's Cilk++ language |17j is a linguistic extension of C++. It is built on the MIT Cilk 
system [5], which is an extension of C. In both of these extensions, programs are ordinary 
programs with a few additional keywords: cilk_spawn, cilk_synk and silk_for. A program 
running on a single processor runs like the original code without the additional keywords. 
The Cilk system contains a work-stealing scheduler. 

• Fastflow [1] is based on lock-free queues explicitly designed for programming streaming 
applications on multi-cores. The authors report that Fastflow exhibits a substantial speedup 
against the state-of-the-art multi-threaded implementation. 

• OpenMP [9] is presented as a set of compiler directives and callable runtime library routines 
that extend Fortran (and separately, C and C++) to express shared memory parallelism. 
It can support pointers and allocatables, and coarse grain parallelism. It also includes a 
callable runtime library with accompanying environment variables. 

• Intel Threading Building Blocks (TBB) |26| is a portable C++ template library for multi- 
core processors. The library simplifies the use of lower level threading packages like Pthreads. 
The library allocates tasks to cores dynamically using a run-time engine. 

• Pthreads [6] is a POSIX standard for the C programming language that defines a large 
number of functions to manage threads. It uses mutual exclusion (mutex) algorithms to 
avoid the simultaneous use of a common resource. 
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2.2. Message Passing Interface systems. There are a few existing libraries that are built on 
top of MPI. These are C++ interfaces for MPI, without job queues, and mostly without higher 
level functionality. It would have been possible to implement MPQueue on top of some of these 
packages, in particular, on top of Boost. MPI. 

• The Boost. MPI library |14[ [To] is a comprehensive and convenient C++ interface to MPI 
with some similarity to our library. Like MPQueue, Boost. MPI uses the BSL to serialize 
messages. 

• MPI++ [2] is one of the earliest C++ interfaces to MPI. The interface is consistent with 
the C interface. 

• mpi++ |15| is another C++ interface to MPI. The interface does not try to be consistent 
with the C++ interface. 

• Object Oriented MPI (OOMPI) [19] is a C++ class library for MPI. It provides MPI func- 
tionality though member functions of objects. 

• Para++ is a generic C++ interface for message passing applications. Like our interface, it 
is a high level interface built on top of MPI and is meant to allow for the quick design of 
parallel applications without a significant drop in performance. The article [8j describes the 
package and includes an example application for PDE. Para++ uses task hierarchies but 
does not implement job queues. 

2.3. Single-System Image (SSI) Cluster software. A single-system image cluster [7j [18] ap- 
pears to be a single system by hiding the distributed nature of resources. Processes can migrate 
between nodes for resource balancing purposes. In some of its implementations, it may be possible 
to run shared memory applications. 

• Kerrighed |20] is an extension of the Linux operating system. It allows for applications using 
OpenMP and Posix multithreading, presenting an alternative to using MPI on a cluster 
with distributed resources. This approach is seemingly simpler than message passing, but 
the memory states of the separate nodes must be synchronized during the execution of a 
shared memory program. This probably is not as efficient as a carefully designed message 
passing solution. 

• OpenSSI is another SSI clustering system based on Linux. It facilitates process migration 
and provides a single process space. OpenSSI can be used to build a robust, high availability 
cluster with no single point of failure. It is a general purpose system that is not specifically 
designed for parallel computing. 

• The Linux Process Migration Infrastructure (LinuxPMI) is a Linux Kernel extension for 
single-system image clustering. The project is a continuation of the abandoned openMosix 
clustering project which is a fork of MOSIX [3]. 

2.4. Summary of Alternatives to MPQueue. As researchers in PDE and combinatorics, we 
have a need for the computational power only obtainable through parallel programming. To be 
practical, we need a method that uses the distributed memory clusters available to us. Our key, 
proven effective programming idea is to use self-submitting job queues. We did not find an exist- 
ing package that completely satisfied all of our requirements, and so wrote our own, MPQueue. 
Our main design goal was to make a simple library that is easy to use, allowing effortless, rapid 
development of scientific experiments requiring parallel execution, without sacrificing performance. 

Existing shared memory model software for multicore systems have some similar functionality 
as MPQueue, but they do not suit all needs because of the limitations on the number of cores on 
currently available hardware. It may be a possibility to use some such systems on top of SSI cluster 
software to implement techniques similar to ours over a distributed memory cluster, but we do not 
believe that this would be any simpler than or have as good performance as our simple, more direct 
message passing approach. 
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FIGURE 3.1. A possible order of jobs in the job queue during the factorization of 
120. The arrows represent the job submission process. The circled jobs are the prime 
number outputs and so they do not produce new submissions. 



Other message passing systems are mainly concerned with lower level tasks. They offer a large 
variety of ways to send and receive messages efficiently but generally do not offer higher level 
constructs. It would probably be a reasonable programming alternative to implement our key idea 
of the self-submitting job queue on top of some of these existing systems, but we do not think 
the result would be as simple to use or have better performance than our implementation. In the 
sequel, we demonstrate with examples how our simple message passing system handles low level 
details automatically and offers an easy to use yet powerful and efficient programming methodology 
based on the self-submitting job queue. 

3. Factoring example 

Listing [3TT1 shows an example program using the MPQueue library for factoring an integer. This 
example demonstrates many of the features of the library. In particular, it shows the general 
structure of a program, including the creation, supervision, and post-processing of job queues. It 
also shows the mechanism by which workers themselves can submit new jobs. The main idea of 
the algorithm is to split the input as a product of two integers that are as large as possible, and 
then submit these factors to the job queue for further splitting. These submissions are done by 
the workers. The job submission process is visualized in Figure 13.11 Note that the order of job 
submissions is not fully determined. 

The companion website |21| contains not only the source code for this and the other examples, 
but a more detailed line- by- line explanation of the code. See also Table [77TI for a complete list of the 
commands found in the MPQueue library, together with a brief description of their function. Here, 
we include an overview of this example program, in particular the lines of code which implement 
key MPQueue features. The MPQueue header file included in all our examples contains our new 
definitions as well as all header files needed to use MPI and the required parts of the BSL |25| . In 
each example, we define a finite list of job types as positive integers, although in this first simple 
example there is only one type of job. 

In each of our example programs, the main function (found here on line 21) is executed by every 
node. First, the MPI is initialized and nodes are split into one boss and several workers. Only the 
boss returns from the MPQstart function call (line 23). The workers are then ready to accept jobs. 
The boss creates two job queues, one for storing jobs to do, and another for storing results. This 
program computes the prime factorization of the example integer found on line 25. The boss places 
the first splitting job into the job queue and starts the supervision of the workers via a call to the 
MPQrunjobs function (line 27), where it will spend the vast majority of its running time. 

In this example, the workers split numbers into factors and submit these factors (see lines 14-15) 
to inqueue for further splitting. The workers return the prime factors, which the boss collects in 
the outqueue. When there are no more jobs, the boss enters a loop (line 29) and retrieves all the 
prime factors from the output queue. 

The user-coded MPQSwitch function is called every time a new job needs to be processed. There 
typically are more than the one job type found in this first most simple example. Here, the workers 
(and only the workers) spend all their productive time in MPQSwitch, recursively finding pairs of 
factors of an input integer and submitting two new factor jobs for each pair of factors found. The 
job. data variable initially contains the serialized input of the job at the time the function is called. 
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#include "MPQueue.h" 
^include "math.h" 
using namespace std; 

const int SPLIT = 1: // only one job type, so no 

// cases in MPQswitch 

void MPQswitch (Tjob k job) { 
int x ; 

from string (x, job. data); 

int sqt = int (sqrt (double (x))); 

for (int y = sqt; y > 1; y ) // search for divisors 

if (0 = x % y) { // found a divisor 

MPQsubmit (Tjob(SPLIT, x/y)); // submit two new jobs 
MPQsubmit (Tjob(SPLIT, y) ) ; 

job . data = ""; // no output to return 

return ; 

} 

} 

int main (int argc , char *argv[]) { 
MPQinit (argc, argv); 

MPQstart (); // only the boss returns 

Tjobqueue inqueue , outqueue; 

int x = 1120581000; // factor this number 

inqueue . push(Tjob (SPLIT , x) ) ; // add one job to the job queue 

MPQrunjobs (inqueue, outqueue); // supervise queue processing 
cout « x« " „ factors w as : \n" ; 

while (! outqueue . empty ()) { // get the results 

int factor; // one factor 

from string( factor , outqueue . front (). data) ; 
cout « factor « " w "; 
outqueue . pop ( ) ; 

} 

MPQstop (); 

} 



LISTING 3.1. Prime factorization example program. 



This variable is replaced by the serialized output of the job. It is a key feature for the ease of doing 
high-level programming with our library that the input and output data can contain arbitrarily 
complicated data structures. In this simple example, the output is the same as the input if no 
factors were found in the loop initiated on line 12, whence the input is a prime, otherwise the 
output is empty. 



4. Matrix square example 

Listing 14.11 shows an example program for calculating the square of a matrix. The example 
demonstrates several new features of the MPQueue library. Namely, it shows how to efficiently 
share a large amount of data with all the workers using MPQsharedata, how to serialize struct data 
types using a template function, and how to use MPQtask to return results to the boss for immediate 
processing rather than using the output job queue. Unlike the first example, this program uses three 
job types so that a switch statement is required in the MPQswitch function. Here, we want the job 
types to be positive integers, so NONE is added to take the unused value of zero. 

The variable matrix contains the input matrix (initialized in lines 43-44), while the variable 
square contains the result of the program, the square of the input matrix. The boss uses our 
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#include "MPQueue.h" 
using namespace std; 

enum { NONE, DATA, MULTIPLY, RESULT }: 
typedef vector < int > Trow; 
vector < Trow > matrix, square; 
typedef struct { 

int pos ; 

Trow result ; 
} Tdata; 



row index 
output row 



template < class Archive > void // needed to serialize Tdata 

serialize (Archive & ar , Tdata & data, const unsigned int version) { 

ar & data . pos ; 

ar & data. result; 

} 

void MPQswitch (Tjob & job) { 
Tdata data; 
switch (job . type) { 
case DATA: 

from string (matrix, job. data); 

break ; 
case MULTIPLY: 

from string (data. pos, job. data); 

data, result = Trow (matrix, size (), 0); 

for (int j = 0: j < matrix, size (); j++) 
for (int k = 0; k < matrix . size (); k++) 

data, result [ j] += matrix [ data . pos ][ k] * matrix [k] [ j ] ; 

job = Tjob(RESULT, data); // prepare the result 

MPQtask (job); // send result to boss 

break ; 

case RESULT: // receive one output row 

from string (data, job. data); 
square [ data . pos ] = data, result; 
job . data = "" ; 

} 



receive the input matrix 

get the row position 
calutate one row 



// nothing to return 



} 



* argv 



{ 







int main (int argc , char 
MPQinit (argc, argv); 
MPQstart (); 

Trow row (10, 1 ) ; 
vector < Trow > mat (row. size 
square . resize (row. size ()); 
MPQsharedata (Tjob(DATA, mat)); 
Tjobqueue inqueue , outqueue; 
for (int i = 0; i < mat. size (); i 

inqueue. push (Tjob (MULTIPLY, i)) 
MPQrunjobs (inqueue, outqueue); 
for (int i = 0; i < row. size (); i 

for (int j = 0; j < row. size () ; 
coufc « square [ i ][ j ] « V; 

cout « "\n" : 

} 

MPQstop (); 

} 



++) 



++) { 
J++) 



input matrix containing 
1 at every entry 
container for the output 
input matrix sent to workers 

every row is a separate job 

run the jobs 

print the output matrix 



LISTING 4.1. Matrix square example program. 
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function MPQsharedata on line 46 to painlessly share an STL structure with the nodes. In lines 48- 
49 the job queue is filled with MULTIPLY jobs, each requesting the calculation of one row of the 
goal matrix. The boss then starts the supervision of the workers, again spending the bulk of its 
time in MPQrunjobs. 

As in the previous example, the workers execute MPQswitch and spend all their productive time 
in that function. Their work, however, this time depends on the job type. Also different from the 
factor example, the MPQswitch function can be executed by the boss, in this case with job. type set 
to RESULT, as the result of a MPQtask call by a worker. The variable data contains one row of 
the result produced by a worker. When the input matrix is shared by all the workers, each worker 
receives the matrix as a DATA job. They deserialize the data and store it locally in the variable 
matrix. Lines 24-29 are where a worker calculates one row of the goal matrix. The input (position 
within the matrix) is deserialized into data.pos, and the calculated row is stored in data. result. 
Workers circumvent the output queue by using MPQtask (line 31) to send their calculated row 
results to the boss. Here, the boss receives and deserializes a row and puts it into the appropriate 
row of the result matrix square. Thus, in lines 33-36, the boss does some actual work, not just 
supervising workers. The boss is usually idle during most of the supervision process. Thus it can 
be more efficient for the boss to do tasks such as post-processing during supervision, rather than 
waiting until the workers are all done. 

5. NON-ATTACKING QUEENS EXAMPLE 

We now present an application that illustrates how to avoid the too many jobs obstacle to scala- 
bility. The example uses local job queues to avoid excessive communication costs. This technique 
enhances the flexibility of our methodology, allowing for the efficient use of the library in somewhat 
unexpected situations. 

The n-queen puzzle is a well-known problem in mathematics. It concerns the placement of n 
non-attacking queens on an n x n chess board. Figure [5Tl shows a valid partial placement with one 
missing queen. For a survey of results on the generalizations of this problem see [4]. Listing [57T1 
shows an example program which counts the number of solutions on an n x n board. Figure 15.21 
shows the serial pseudo code. We were able to run this code for n < 20 on a cluster containing 24 
Intel(R) Xeon(TM) 2.80GHz dual CPUs with hyper-threading. For n = 20, it took 9.0 hours (see 
Figure E3|) using the 96 = 24 x 2 x 2 cores to find the 39, 029, 188, 884 solutions. Note that the 
state of the art is n = 26, that is, the number of solutions for n = 27 is not known. 

Our code uses a parallel version of a standard backtracking algorithm. A worker keeps a local job 
queue containing possible search branches, i.e., valid partial placements. If this queue grows large 
enough, then the worker submits one of the partial placements to the boss node as a job so that 
another worker can process it. Submitting jobs to the boss node too early may result in a very large 
global job queue and a lot of unnecessary communication between the nodes. Submitting jobs too 
late can starve the workers for jobs. The decision depends on the number of nodes, the speed of the 
nodes, and the branching factor of the search tree. Our code uses a simple heuristic that depends 
on an overflow value that we have adjusted experimentally. Figures [5.3[ I5.4[ and 15.51 show the effect 
of different choices of the value. This very simple submission process could be fine tuned using the 
MPQinfo command (see Table [7TT]) . the only MPQueue command not demonstrated in any of our 
three examples. Many of the task distribution schemes discussed in the survey paper |13| could be 
implemented using the MPQinfo command. 

In the sequel of this section we provide a description of the more important aspects of the 
algorithm as they are specifically implemented in the lines of code found in Listing 15.11 The size 
of the board is set to some positive integer, 20 in this example (line 9). A partial placement of 
queens is a vector of type Trow of length between and size, whose integer value j in the i-th. entry 
denotes a valid (non-conflicting) placement of a queen in the j-th row of the i-th column of the 
board. When the length of a partial placement row equals size, it is in fact a full placement and a 
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#include "MPQueue.h" 
^include <deque> 
using namespace std; 

enum { NONE, PLACE, RESULT }: 
typedef vector < int > Trow; 
unsigned long int allsolutions 
int overflow = 70; 
const int size = 20; 



0; 



total number solutions 
contr oils local queue size 
size of the board 



bool inline fits (const Trow & row) { 
int j = row. size () — 1; 
for (int i = 0; i < j; 

if ( ( i'ow[i] = row. back () ) 
return false ; 
return true ; 

} 



job) { 



( abs 



- row\j } ) = j - 
queens interfere 
last queen fits 



i ) ) 



0: 



// a partial placement 

// local job queue 

// local number of solutions 

// add queen in next column 

// populate a local job queue 

// still local jobs to do 



void MPQswitch (Tjob 
Trow row; 

deque < Trow > rows; 
unsigned int solutions 
switch (job. type) { 
case PLACE: 

from string (row, job. data); 
rows . push_back (row); 
while (! rows . empty ()) { 
row = rows . back ( ) ; 
rows . pop _back (); 

for (row . push _ back (0); row. back () < size; row. back () + + ) 

if (fits (row)) // does the new queen fit? 

if (row. size () = size) // all queens added 

solutions^ — h; // found a new solution 

else { 

rows . push_back (row); // add to local job queue 

if (rows, size () > overflow) { // if too many local jobs 

MPQsubmit (Tjob(PLACE, rows, front ())); // send one to the boss 
rows . pop_ front (); 



} 



} 



} 



} 



job = Tjob (RESULT, solutions); 
MPQtask (job); 
break ; 
case RESULT: 

from string (solutions, job . data) 
allsolutions += solutions ; 
job . data = "" ; 

} 



// prepare the result 

// send result to the boss 



// update total solutions 

// with new result 

// no result to return 



int main (int argc , char *argv[]) { 
MPQinit (argc, argv); 
MPQstart (); 

Tjobqueue inqueue , outqueue; 

inqueue . push (Tjob (PLACE, Trow ())); 

MPQrunjobs (inqueue, outqueue); 

cout « allsolutions « " w solutions\n" 

MPQstop (); 

} 



// initial empty placement 



LISTING 5.1. Non-attacking queens example program. 
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FIGURE 5.1. A partial placement of queens on a 5 x 5 board. In the example code, 
this placement is encoded as the vector (0,3,1,4). This partial placement can be 
extended to the full placement (0, 3, 1, 4, 2). There are 10 distinct placements on the 
5x5 board. 



Input: size 
Output: solutions 



add empty row to the queue of partial placements 
while queue is not empty do 

move the last job from queue into row 
increase the size of row 
for i E {0, . . . , size} do 

set the last coordinate of row to i / 
if row is a good partial placement then 
if row is a full placement then 
increment solutions 

else 

add row to queue 



no queens placed yet 



add room for the next column 



place the next queen into the z-th row 



FIGURE 5.2. Pseudo code for the serial non-attacking queens placement program. 
The algorithm uses simple back tracking. 



solution has been found. After the usual initialization where the boss and worker nodes are created, 
the boss puts an empty placement job on the queue and begins supervision by entering MPQrunjobs 
(line 57), where it will stay until the counting is done and its global variable allsolutions is output 
with the final result. In the algorithm's counting of numbers of solutions at the ends of branches 
of partial placements, workers return partial sums counting some valid placements to the boss by 
submitting a RESULTS task to the boss via a call to MPQtask (line 43). This causes the boss to 
update its global variable of the master count (line 47). 

In order to control communication costs (and in fact fine-tune performance), workers use local 
job queues named rows (line 27) containing up to overflow (line 8) number of partial placement row 
vectors. Until the queue is empty, workers pop off local partial placement jobs, repeatedly calling 
the function fits (lines 11 and 31) with a row vector containing an attempt to extend to a valid 
placement to the next column. A valid extension causes the worker to either increment the number 
of solutions found (line 33) if the edge of the board was reached, or otherwise add it as a partial 
placement job to the local queue for another round of extending by itself (line 35). Before the next 
local job is started, the worker checks if the local queue has reached overflow, and if so pops a partial 
placement job off of the local queue and sends it to the boss to be placed on the global queue (lines 
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overflow = 70 overflow = 85 




time (sec) time (sec) 



FIGURE 5.3. Load diagrams of the 20-queen placement example code for two dif- 
ferent values of the overflow variable. In both cases there were 95 available workers. 
The shaded region shows the number of workers actively working on jobs. The solid 
curve shows the number of jobs available in the global job queue together with the 
number of active workers. The unshaded region below 95 represents idle workers. 
The pictures show that the choice of 85 for overflow is too large, since with this 
choice there are not enough available jobs and so the program takes longer to finish. 
The optimal value of 70 was determined experimentally. From Figures [5.41 and T5.5I it 
appears that the optimal overflow value is somewhat predictable and depends only 
on the size of the problem, not the size of the cluster. 



Nodes p 


Run time T(p) 


Worker usage 


Total usage 


Efficiency p ^ ( ] p) 


1 


94.08 sec 




94.08 CPU sec 




3 


47.05 sec 


94.10 CPU sec 


141.15 CPU sec 


0.67 


6 


18.77 sec 


93.85 CPU sec 


112.62 CPU sec 


0.84 


24 


4.11 sec 


94.53 CPU sec 


98.64 CPU sec 


0.95 


48 


2.05 sec 


96.35 CPU sec 


98.40 CPU sec 


0.96 


72 


1.44 sec 


102.24 CPU sec 


103.68 CPU sec 


0.91 


96 


1.11 sec 


105.45 CPU sec 


106.56 CPU sec 


0.88 



TABLE 5.1. Processing times for the 15-queen placement problem using the overflow 
value 30. The Worker usage column is the number of workers multiplied by the run 
time, while Total usage includes the boss as well. A constant worker usage would 
indicate perfect scaling. For larger board sizes, the efficiency increases up to and 
beyond our maximum available 96 nodes. 



36-38). When a worker's local queue is finally empty, the worker returns the partial result via a 
call to MPQtask and then terminates (lines 43-44). The worker is of course then available for the 
boss to give it another partial placement job off of the global queue, until that queue is itself empty 
and the final count has been compiled. 

Table [57T1 and Figures [5.41 and [531 show that our approach to the problem scales well. For the size 
15 board, we could have used more than the available 96 nodes; to attempt the unknown n = 27 
case we could make good use of a very large number of nodes. 
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40 50 60 
nodes p 

FIGURE 5.4. Speedup as a function of nodes for the 12-queen placement problem. 
The bold numbers indicate the overflow values used for a given curve. The dashed 
line is the theoretical maximum speedup of p — 1, the number of workers. It appears 
that adding more than 55 nodes does not increase the speedup. This limitation is 
the result of the small board size. The number of solutions is only 14,200. 




nodes p 



FIGURE 5.5. Speedup as a function of nodes for the 15-queen placement problem. 
The bold numbers indicate the overflow values used for a given curve. It appears 
that the optimal overflow value does not depend on the number of nodes, only on the 
size of the problem. The dashed line again corresponds to the theoretical maximum 
speedup. The number of solutions for this board size is 2,279,184. 
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We did not take advantage of the symmetry of the problem, which would have immediately 
resulted in a four-fold speedup. This could be done with minimal work, but our goal here is to 
present a meaningful yet simple example demonstrating the use of our library and methodology. A 
possible future extension to the MPQueue library that could help reduce the number of idle workers 
for problems like this one would be the implementation of priority job queues. 

As a reminder, the companion website |21j contains the source codes and corresponding detailed 
line-by-line explanations for this and the previous two examples, along with the MPQueue library 
and the few other files needed to compile and execute the programs. Table [7+1 contains a complete 
list of MPQueue commands, along with a brief description of their function. 

6. Application to Nonlinear Elliptic PDE 

Our library has provided us with an easy way to port to a parallel environment our serial code 
for solving PDE. The resulting increase in computational power has enabled us to solve the com- 
putationally intensive problem 

(6.1) Au + f s (u) = inn 

u = on d£l, 

where A is the Laplacian, is the square or cube, and f s is the family of nonlinearities f s (u) = 
su + u 3 , parameterized by s € R. Other regions and nonlinearities can be handled as well. We 
had previously developed serial C++ code for obtaining good approximations of solutions to pa- 
rameterized PDE such as (I6.ip . provided that the dimensions involved were not too big. The reader 
can refer to |22[ [23] (serial) and |24| (parallel) for the mathematical details of our PDE algorithms, 
which are based on Newton's method operating in a coefficient space corresponding to eigenfunction 
expansions of solution approximations. Generally, our C++ code follows branches of solutions by 
starting with an initial point on a branch and an approximate tangent vector to the branch, applying 
Newton's method to find a next point and corresponding tangent, and repeating until a window is 
exited. Each time a new bifurcation point is identified on a branch, our code seeks points on new, 
bifurcating branches, which can then also be followed. Processing a bifurcation point can be very 
quick but it can also take more time than following a branch. Thus, it is natural that following a 
branch and analyzing a bifurcation point are the two job types we use. 

The required input for both job types is complicated. Each input is a structure with 18 dif- 
ferent fields, containing such data as the discretized solution approximation at N grid points, the 
corresponding M eigenfunction expansion coefficients, the parameter s, the current tangent vector, 
the branch history, and C++ parameter values that lead to the current state (speed, tolerances, 
maximum iteration counts, etc.). For large problems, the embedded vectors may be very large. The 
automatic serialization feature of our library makes it easy to pass such data between nodes. 

Figure [6~T1 shows a partial bifurcation diagram when £1 is the square (0, ir) 2 . The horizontal axis is 
the parameter s and the vertical axis is the value of the solution at a generic point (x* , y*) |23| . The 
diagram demonstrates how branch following creates a tree structure whose growth is unpredictable. 
Details of the size and speed for this simple example are shown in the first row of Table 16.11 Using 
two nodes is essentially a simulation of our serial code, since there is only one worker. 

The second row of Table [6+1 shows the same summary for a problem where f2 is the cube (0, 7r) 3 
and we search for the branches connected to the bifurcation point on the trivial branch at s = 
18. Parallelization is essential since it takes about a minute to compute each of the thousands of 
solutions. 

We conclude the section with a few details of our implementation for creating bifurcation di- 
agrams. Solution approximations lie in a subspace spanned by the first M eigenfunctions of the 
Laplacian, which themselves have been previously and independently approximated by M vectors 
in R , obtained via calls to ARPACK if not known in closed form. For the example results included 
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1.5 



0.5 



-0.5 



ft: 




b-2 



Pa 



P 5 



FIGURE 6.1. The first four primary branches and all the branches connected to them 
for PDE (|6.ip on a square region. For this example, a total of 12 jobs are placed 
on the job queue. The lower diagram represents a possible creation order of these 
jobs in the job queue during the generation of the bifurcation diagram presented in 
the upper diagram. A solid arrow represents the submission of a bifurcation analysis 
job, while a dotted arrow represents a branch following job. The trivial branch u = 
lying on the horizontal axis is the first branch to be followed. This job, labeled b±, 
is the first and only job to be placed in the job queue by the boss. One of the 
workers follows this trivial branch and encounters the 3 primary bifurcation points 
and submits 3 bifurcation jobs P\,P2, and P3. These three bifurcation points are 
each processed by a different worker, which leads to the submission of 4 branch jobs. 
This process continues until the job queue is empty. The workers do not send results 
to the boss using the outqueue; they store all their results in separate files which are 
post processed by other scripts after the termination of the program. 



in this section, these bases are well known in terms of sine functions. The construction of the Ja- 
cobian of the object function for Newton's method requires order M 2 numerical integrations, each 
requiring order N arithmetic operations. The number of integrations is reduced somewhat in the 
presence of symmetry. The search direction system is solved via a standard LAPACK subroutine. 
Each approximated point requires roughly 4 iterations of Newton's method. Finding bifurcation 
points requires the computation of the eigenvalues of the Jacobian, computed by another LAPACK 
routine. 

7. IMPLEMENTATION 

7.1. Design. Passing complicated data structures between nodes requires a significant amount of 
work in MPI. To avoid this difficulty, we rely on the BSL. This library encodes every standard STL 
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o 


M 


N 


Branches 


Bifurcations 


Solutions 


Nodes 


Run time 


Worker usage 


(0,vr) 2 


323 


41 2 


7 


5 


165 


2 
4 


8.17 min 
3.73 min 


8.17 CPU min 
11.19 CPU min 




564 


21 3 


385 


364 


19724 


26 
50 
76 


22.4 hour 

11.5 hour 
8.98 hour 


560.0 CPU hour 
563.5 CPU hour 
673.5 CPU hour 



TABLE 6.1. Processing times for PDE examples. These results are using NCSA's 
IA-64 TeraGrid Linux Cluster (Mercury) with 1.3 GHz nodes |11] , The data for 
the square region corresponds to the diagrams found in Figure 16.11 The entries for 
the cube region are typical of the results found in [24]. Note that the cube problem 
scales well since there are many jobs. 



data structure automatically as a string, and more complicated data types can be encoded with a 
minimal amount of work. The serialization is done mostly automatically using template functions. 
A serialized data string is sent between the nodes in two steps. In the first step, a preliminary 
message is sent containing the size of the string together with some additional information like the 
job type. In the second step, the string itself is sent. Most of this is done using point-to-point 
communication with MPI_Send and MPI_Recv. The exception is the MPQsharedata command, 
which uses point-to-point communication for the preliminary message and uses M PI _ Beast for the 
data string. 



Input: none 
Output: none 

repeat 

wait for a message from the boss 
switch the message is a 
case data share request 

accept the broadcast message from the boss 
call MPQswitch to handle the broadcast message 

case job 

call MPQswitch to run the job 
send the result back to the boss 

case stop command 
stop 



FIGURE 7.1. Pseudo code for the main loop of the workers. The workers get to 
this loop from the MPQstart function. 

As a result of calling MPQstart, all the nodes become workers except node zero, which becomes 
the boss. The workers go into a waiting loop, as shown in Figure 17.11 and the boss continues its 
own work. When the boss requires help, it builds a queue of jobs and becomes a supervisor using 
MPQrunjobs, as shown in Figure 17.21 The boss then assigns jobs in the job queue to available 
workers and removes these jobs from the job queue. The workers execute MPQswitch with the job. 

When a worker finishes a job, it sends a message to the boss, letting the boss know that the 
result is available. The boss accepts the result from the worker and stores it in the output queue. 
The boss keeps track of the number of jobs sent out. If the job queue is empty or every worker is 
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Input: input queue 
Output: output queue 



while input queue is not empty or workers are working do 

if input queue is not empty and there are idle workers then 
assign a job to a worker 
remove the job from the input queue 
else 

wait for a message from a worker 
switch the message is a 
case job submission 

add the job to the input queue 

case job result 

if result is not empty then 

add the result to the output queue 

case task 

call MPQswitch to run the task 
send the result back to the worker 

case info request 

send the info back to the worker 



Figure 7.2. Pseudo code for the MPQrunjobs function. This function is run by the boss. 



working already, then it only accepts messages and does not try to assign jobs. The supervision 
phase ends when the job queue is empty and the result of every assigned job is received. 

During the supervision phase, the workers can create new jobs by sending a job to the boss using 
MPQsubmit. The boss puts these jobs into the currently executing job queue. 

The workers can use MPQtask to ask the boss to execute a task immediately and return the result 
to the same worker. To satisfy this request, the boss executes MPQswitch. This feature can be used, 
for example, to set or get the value of a global variable like a counter. This allows for rudimentary 
communication between the workers. The feature can also be used for returning results to the 
boss. This allows the boss to post-process the results while the workers are busy with their jobs. A 
sequence diagram of the supervision phase is shown in Figure 17.31 At the end of the program, the 
boss uses MPQstop to tell the workers to quit. 

See Table 17.11 for a complete list of the commands found in the MPQueue library, together with 
a brief description of their function. The interested reader can consult the companion website |21| 
to find a detailed reference guide. 

7.2. Scaling and overhead. In the non-attacking queens example we presented evidence suggest- 
ing that our methodology scales well. In particular, we show with load diagrams that all workers 
can be kept effectively busy if the problem has very many small jobs that can be managed via a local 
job queue to avoid excessive communication. We show that the efficiency of our parallel solution 
to this typical application is very good, up to a certain number of nodes, and that this number of 
nodes is very large if the problem is sufficiently big. We show that the speedup for this example 
can be close to ideal when using the optimal overflow, that this optimal overflow value appears to 
be independent of the number of nodes and only dependent on the size of the problem, and that 
for a sufficiently large problem, an arbitrary number of nodes can be effectively used. For our PDE 
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nodeO:boss 



node 1: worker 



node2: worker 



'data (MPQsharedata) 



MPQswitch is running 1 ^ * 



^ I job queue (MPQrunjobs) 
job 



(MPQinfo) 



info 



job 



jobqueue item (MPQsubmit) 



if 



MPQswitch is running 1 ^ 



task result 



task (MPQtask) 



MPQswitch is running 1 ^ 



task result 



job result 



MPQswitch is running 1 ^ . 



task (MPQtask) 



_ job result queue 



job result 



FIGURE 7.3. An example sequence diagram for the communication between the 
boss and worker nodes. The boss sends out some global data to the workers using 
MPQsharedata. The boss creates a job queue with a single job and then becomes 
the supervisor by calling MPQrunjobs. The first worker takes this job and realizes 
that an extra job should be created. It uses MPQinfo to check if there are available 
workers to take this additional job. The second worker is not busy, so the first worker 
submits the new job. The new job is assigned to the second worker. Both workers 
need the value of a global variable so they use MPQtask to get it from the boss. 
When the workers are finished, they send the results to the boss and the supervision 
phase stops. 
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struct Tjob {int type; string data; }; 


Define job data type. 


typeder queue < ljob> ljobqueue; 


Denne jobqueue data type. 


void MPQinit(int argc, char *argv[ ]); 


Initialize MPI. 


void MPQstart(); 


Split the workers from the boss process. 


void inline MPQsubmit(const Tjob & job): 


Submit a job into the currently running job queue. 


template <class T> void 
inline MPOtaskfT'ob &- iobV 


Ask the boss to run a task. 


void MPQinfo 

(int & queuesize, int & sent); 


Get the number of jobs in the running job queue 
and the number of workers currently working. 


void MPQmnjobs 

(Tjobqueue & inq, Tjobqueue & oufcq); 


Assign the jobs in inq to the workers and collect 
the results in outq. 


void inline MPQsharedata(const Tjob & job); 


Send data to all the workers. 


void MPQswitch( Tjob & job); 


Execute a job. 


void MPQstop(); 


Release the workers and stop MPI gracefully. 


template < class T> string 
to string(const T&in); 


Serialize any variable. 


template <class T> void 

from string(T & out, const string & str); 


Deserialize a variable. 


#define LOAD_DIAGRAM 


Generate load data. 



Table 7.1. A complete list of the commands found in the MPQueue library. 



application where there are a lesser number of much larger jobs and local job queues are not needed, 
we have observed similarly good performance indicators (see Table [6TT1 and |24|). 

In a final test, we demonstrate that the overhead of our implementation is reasonable. To observe 
the cost of bookkeeping, data serialization, and communication, we ran two experiments where each 
job was an instruction to the workers to sleep for one second. In the first experiment, we sent no 
additional input or output data, only the type of the job. In the second experiment, each job was 
sent with a vector containing 1000 doubles as input, and then the same vector was returned as 
output. The addition of data causes some unavoidable delay on our distributed memory system 
due to the serialization time and inherent communication speed of the network. The results are 
summarized in Table 17.21 It shows that the correspondence between the overhead and the number 
of jobs is linear in both experiments. 

8. Conclusions 

We have demonstrated that parallel self-submitting job queues are a powerful computational 
paradigm. Job queues are high-level constructs not available in other MPI libraries. The ability of 
workers to submit jobs to the job queue is the key component that gives our method its flexibility 
and power. We implemented our ideas in a light-weight and easy to use library based on MPI. The 
library, MPQueue, has been used successfully in real research applications. Using our interface saves 
effort by avoiding low-level coding; scientists without expertise in parallel programming can rapidly 
develop code to obtain good results for serious research problems. Experts in parallel programming 
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without data 


with data 


Jobs n 


Run time T(p) 


Overhead T(p) — n/95 


Run time T(p) 


Overhead T(p) — n/95 


760 


8.02 


0.02 


8.73 


0.73 


3040 


32.07 


0.07 


24.88 


2.88 


12160 


128.19 


0.19 


139.50 


11.50 


48640 


512.76 


0.76 


557.90 


45.90 


194560 


2051.09 


3.09 


2231.28 


183.28 


778240 


8204.51 


12.51 


8964.35 


772.35 



TABLE 7.2. The time in seconds to execute sleep(l) jobs on 96 nodes. We used the 
same cluster specified in Section [5l Since there are 95 workers, it takes at least n/95 
seconds to do n such jobs. The rest of the time is the overhead. The overhead is 
approximately 16 microseconds per job with no data sent, and 98 milliseconds per 
job with a vector of 1000 doubles serialized, deserialized and sent both ways for every 
job. 



can also take advantage of our methodology for applications that can be broken into relatively large 
pieces. In applications such as our PDE example found in Section [6j this decomposition is obvious. 
We also find the approach to be effective in less obvious situations, as demonstrated in the queens 
example in Section [5] which uses local queues on top of the global job queue to avoid the excessive 
communication that would result from transmitting too many small jobs across nodes. We have 
supplied evidence in this same example that our approach is scalable. The overhead associated with 
our implementation has been demonstrated to be reasonable. Possibilities for applications are wide; 
we have also used MPQueue to find winning strategies in combinatorial game theory |27| and in 
simulating a critical branching random walk [1Q]. 

Further improvements to MPQueue could include developing algorithms for automatically ad- 
justing the local queue size (i.e., overflow), or other more sophisticated control procedures for dis- 
tributing workload, such as the implementation of priority job queues and the capability of workers 
to send jobs directly to workers. We could easily add to the library the ability for a subset of work- 
ers to act as sub-supervisors, each with their own pool of workers. We believe that the currently 
implemented features of MPQueue will suffice for most applications, and that its simplicity is in 
many ways an advantage. 

In some communication intensive applications, efficiency could be increased by avoiding the se- 
rialization of fixed- length data. The Boost. MPI library sends fixed length data more efficiently, 
without sacrificing convenience. It would be possible to implement MPQueue on top of Boost. MPI 
or other MPI interfaces in order to take advantage of the optimizations offered by those libraries. 
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